program cubeflu

use basic_var
use mpi_basic_var
use timing
use mpi
#ifdef _ACCEL
use accel_lib
#endif
#ifdef _OPENMP
use omp_lib
#endif

real(myk) :: start_time,end_time,elapsed_time
real(myk) :: step_start_time,step_end_time,step_elapsed_time=0._myk
real(myk) :: io_step_start_time,io_step_end_time,io_step_elapsed_time=0._myk
integer :: ik

call start_mpi()

GPU = .false.

#ifdef _ACCEL
print*,'NEW n_rank,NCOREXNODE: ',n_rank,NCOREXNODE
if(mod(n_rank,NCOREXNODE) .lt. NGPUXNODE) then
  call acc_set_device(acc_device_nvidia)
  device_kind = acc_get_device()
  print*,'Selected device kind:  ',device_kind
  num_devices = acc_get_num_devices(acc_device_nvidia)
  print*,'Number of devices available: ',num_devices
  call acc_set_device_num(mod(n_rank,NCOREXNODE),acc_device_nvidia)
  print*,'n_rank: ',n_rank,' tries to set GPU: ',mod(n_rank,NCOREXNODE)
  my_device = acc_get_device_num(acc_device_nvidia)
  print*,'n_rank: ',n_rank,' is using device: ',my_device
  print*,'Set GPU to true for rank: ',n_rank
  GPU = .true.
endif
#endif
#ifdef _OPENMP
if(.not. GPU) then
   if(openmp_activate == 1) then
     print*,'Setting OpenMP threads :',NCOREXNODE - NGPUXNODE
     call omp_set_num_threads(NCOREXNODE - NGPUXNODE) ; NCOREXNODE = NGPUXNODE + 1
   else ; call omp_set_num_threads(1) ; endif
endif
#else
if(openmp_activate == 1) then
  print*,'Error: please compile using OpenMP or disable openmp_actiavte in mpi_gpu.dat' ; call end_mpi() ; STOP
endif
#endif

call time_init()

call size_init_001()
call decompose_mpi_001()
call allocate_storage_001()
call random_init_001()
call input_grid_001()
call input_grid_tot_001()
call ini_eval_001()

call size_init_002()
call decompose_mpi_002()
call allocate_storage_002()
call random_init_002()
call input_grid_002()
call input_grid_tot_002()
call ini_eval_002()

call allocate_storage_exchange_001()
call allocate_storage_exchange_002()
print*,'MPI and Memory allocated succesfully'

call read_scalars_001()
call read_scalars_002()

if(itmin.eq.0) then
   call var_init_001()
   call var_init_002()
   !call print_sol_binary_001() ! da testare
   !call print_sol_binary_002() ! da testare
   call write_sol_binary_pio_001()
   call write_sol_binary_pio_002()
   call print_vtk_001()
   call print_vtk_002()
else   ! restart
   !call read_sol_binary_001
   !call read_sol_binary_002
   call read_sol_binary_pio_001()
   call read_sol_binary_pio_002()
   t = itmin*dt
endif

#ifdef POST_RUNTIME
if(itmin.lt.post_itmin) then
   call init_pos_001()
   call init_pos_002()
else
   call read_average_tecplot_space_and_time_001(1)
   call read_average_tecplot_space_and_time_002(1)
endif
#endif

call var_k_init_001()
call rhs_k_init_001()

call var_k_init_002()
call rhs_k_init_002()

if(GPU) then
  call update_gpu_001()
  call update_gpu_002()
  print*,'GPU activated for rank : ',n_rank
endif
call print_screen_001()
call print_screen_002()

call mpi_barrier(mpi_comm_world,ierr)
call time(start_time)
do it=itmin+1,itmax
   call time(step_start_time)
   t = t+dt
   if(n_rank.eq.0) then
     if(mod(it,itout).eq.0) then
       print*, '__________________________________________________'
       print*, 'it = ',it,' t = ',t
     endif
   endif
   !azzeramento rhs per il prossimo passo
   if(GPU) then ; call rhs_init_acc_001() ;       else ; call rhs_init_omp_001() ;       endif
   if(GPU) then ; call rhs_init_acc_002() ;       else ; call rhs_init_omp_002() ;       endif
   do ik=1,s
!      print*,'debug ik: ',ik
      tk = t + alpha(ik)*dt
      ! First block: step I
      if(GPU) then ; call var_k_eval_acc_001(ik) ;       else ; call var_k_eval_omp_001(ik) ;       endif
      if(GPU) then ; call update_var_k_mpi_acc_001() ;   else ; call update_var_k_mpi_omp_001() ;   endif
      if(GPU) then ; call exc_bc_var_k_acc_001();        else ; call exc_bc_var_k_omp_001() ;       endif

      ! Second block: step I
      if(GPU) then ; call var_k_eval_acc_002(ik) ;       else ; call var_k_eval_omp_002(ik) ;       endif
      if(GPU) then ; call update_var_k_mpi_acc_002() ;   else ; call update_var_k_mpi_omp_002() ;   endif
      if(GPU) then ; call exc_bc_var_k_acc_002();        else ; call exc_bc_var_k_omp_002() ;       endif

      ! Wait exchange var
      if(GPU) then ; call exc_wait_var_k_acc_001();      else ; call exc_wait_var_k_omp_001() ;       endif
      if(GPU) then ; call exc_wait_var_k_acc_002();      else ; call exc_wait_var_k_omp_002() ;       endif

      ! First block: step II
      if(GPU) then ; call bc_var_k_acc_001() ;           else ; call bc_var_k_omp_001() ;           endif
      if(GPU) then ; call rhs_k_eval_acc_001() ;         else ; call rhs_k_eval_omp_001() ;         endif
      if(GPU) then ; call exc_bc_rhs_k_acc_001();        else ; call exc_bc_rhs_k_omp_001() ;       endif

      ! Second block: step II
      if(GPU) then ; call bc_var_k_acc_002() ;           else ; call bc_var_k_omp_002() ;           endif
      if(GPU) then ; call rhs_k_eval_acc_002() ;         else ; call rhs_k_eval_omp_002() ;         endif
      if(GPU) then ; call exc_bc_rhs_k_acc_002();        else ; call exc_bc_rhs_k_omp_002() ;       endif

      ! Wait exchange rhs
      if(GPU) then ; call exc_wait_rhs_k_acc_001();      else ; call exc_wait_rhs_k_omp_001() ;       endif
      if(GPU) then ; call exc_wait_rhs_k_acc_002();      else ; call exc_wait_rhs_k_omp_002() ;       endif

      ! First block: step III
      if(GPU) then ; call bc_rhs_k_acc_001() ;           else ; call bc_rhs_k_omp_001() ;           endif
      if(GPU) then ; call rhs_eval_acc_001(ik) ;         else ; call rhs_eval_omp_001(ik) ;         endif

      ! First block: step III
      if(GPU) then ; call bc_rhs_k_acc_002() ;           else ; call bc_rhs_k_omp_002() ;           endif
      if(GPU) then ; call rhs_eval_acc_002(ik) ;         else ; call rhs_eval_omp_002(ik) ;         endif
   enddo
   ! First block
!HARMFUL   if(GPU) then ; call bc_rhs_acc_001() ;            else ; call bc_rhs_omp_001() ;               endif
   if(GPU) then ; call var_eval_acc_001() ;          else ; call var_eval_omp_001() ;             endif
   if(GPU) then ; call update_var_mpi_acc_001() ;    else ; call update_var_mpi_omp_001() ;  endif
!   if(GPU) then ; call bc_var_acc_001() ;            else ; call bc_var_001_omp() ;          endif ! DOPPIONE PERCHE' FATTA AL I SUB-STEP RK
!   if(GPU) then ; call get_exc_bc_acc_001();         else ; call get_exc_bc_001_omp()   ;    endif ! DOPPIONE PERCHE' FATTA AL I SUB-STEP RK
!   if(GPU) then ; call put_exc_bc_acc_002on001() ;   else ; call put_exc_bc_002on001_omp() ; endif ! DOPPIONE PERCHE' FATTA AL I SUB-STEP RK
   if(GPU) then ; call post_eval_acc_001() ;         else ; call post_eval_omp_001() ;       endif

   ! Second block
!HARMFUL   if(GPU) then ; call bc_rhs_acc_002() ;            else ; call bc_rhs_omp_002() ;          endif
   if(GPU) then ; call var_eval_acc_002() ;          else ; call var_eval_omp_002() ;        endif
   if(GPU) then ; call update_var_mpi_acc_002() ;    else ; call update_var_mpi_omp_002() ;  endif
!   if(GPU) then ; call bc_var_acc_002() ;            else ; call bc_var_002_omp() ;          endif ! DOPPIONE PERCHE' FATTA AL I SUB-STEP RK
!   if(GPU) then ; call get_exc_bc_acc_002();         else ; call get_exc_bc_acc_002_omp() ;  endif ! DOPPIONE PERCHE' FATTA AL I SUB-STEP RK
!   if(GPU) then ; call put_exc_bc_acc_001on002() ;   else ; call put_exc_bc_001on002_omp() ; endif ! DOPPIONE PERCHE' FATTA AL I SUB-STEP RK

   if(GPU) then ; call post_eval_acc_002() ;         else ; call post_eval_omp_002() ;       endif

   call time(step_end_time)
   step_elapsed_time = step_elapsed_time + (step_end_time - step_start_time)
   if(mod(it,itout).eq.0) then
     call time(io_step_start_time)
     if(GPU) call update_cpu_001()
     if(GPU) call update_cpu_002()
     call print_screen_001()
     call print_screen_002()
!     call print_sol()
!     call print_sol_binary_001()
!     call print_sol_binary_002()
     call write_sol_binary_pio_001()
     call write_sol_binary_pio_002()
     call print_vtk_001()
     call print_vtk_002()
     call time(io_step_end_time)
     io_step_elapsed_time = io_step_end_time - io_step_start_time
     if(n_rank.eq.0) then
       open(unit=1,file='iter.dat')
         write(1,*) it
         write(1,*) it + itmax-itmin
         write(1,*) itout
         write(1,*) dt
#ifdef POST_RUNTIME
         write(1,*) post_itmin  ;  write(1,*) it + itmax-itmin  ;  write(1,*) post_itout
#endif
       close(1)
     endif
     if(n_rank.eq.0) print*,'Step Elapsed time : ',step_elapsed_time/itout ; step_elapsed_time = 0._myk
     if(n_rank.eq.0) print*,'IO/Step Elapsed time : ',io_step_elapsed_time 
   endif
#ifdef POST_RUNTIME
   post_it = it
   if((it .ge. post_itmin) .and. (mod(it,post_itout).eq.0)) then
     call time(post_step_start_time)
     if(GPU) call update_cpu_001()
     if(GPU) call update_cpu_002()
     call print_screen_001()                      ;  call print_screen_002()   

     call enlarge_mpi_space_and_time_first_001()  ;  call enlarge_mpi_space_and_time_first_002()  
     call post_mpi_first_001()                    ;  call post_mpi_first_002()  
     call reduce_mpi_space_and_time_first_001()   ;  call reduce_mpi_space_and_time_first_002()  

     call enlarge_mpi_space_and_time_second_001() ;  call enlarge_mpi_space_and_time_second_002()  
     call init_pos_space_average_001()            ;  call init_pos_space_average_002()  
     call post_mpi_second_001()                   ;  call post_mpi_second_002()  
     call reduce_mpi_space_second_001()           ;  call reduce_mpi_space_second_002()  
     call reduce_mpi_space_and_time_second_001()  ;  call reduce_mpi_space_and_time_second_002()  
     if(n_rank .eq. 0) call print_average_tecplot_space_001(1)
     if(n_rank .eq. 0) call print_average_tecplot_space_002(1)

     call time(post_step_end_time)
     post_step_elapsed_time = post_step_end_time - post_step_start_time
     if(n_rank.eq.0) print*,'POST/Step Elapsed time : ',post_step_elapsed_time 
   endif
#endif
enddo

it = itmax ! otherwise it=itmax+1 
!call print_sol_binary()
#ifdef POST_RUNTIME
if(n_rank .eq. 0) call print_average_tecplot_space_and_time_001(1)
if(n_rank .eq. 0) call print_average_tecplot_space_and_time_002(1)
#endif

call time(end_time)
!print*,'ending Runge-Kutta'
elapsed_time = end_time - start_time
print*,'Elapsed time : ',elapsed_time,' GPU? ',GPU

call end_mpi()

end program cubeflu
